QTL mapping for resistance against cereal cyst nematode (Heterodera avenae Woll.) in wheat (Triticum aestivum L.)

The resistance to cereal cyst nematode (Heterodera avenae Woll.) in wheat (Triticum aestivum L.) was studied using 114 doubled haploid lines from a novel ITMI mapping population. These lines were screened for nematode infestation in a controlled environment for two years. QTL-mapping analyses were performed across two years (Y1 and Y2) as well as combining two years (CY) data. On the 114 lines that were screened, a total of 2,736 data points (genotype, batch or years, and replication combinations) were acquired. For QTL analysis, 12,093 markers (11,678 SNPs and 415 SSRs markers) were used, after filtering the genotypic data, for the QTL mapping. Composite interval mapping, using Haley-Knott regression (hk) method in R/QTL, was used for QTL analysis. In total, 19 QTLs were detected out of which 13 were novel and six were found to be colocalized or nearby to previously reported Cre genes, QTLs or MTAs for H. avenae or H. filipjevi. Nine QTLs were detected across all three groups (Y1, Y2 and CY) including a significant QTL "QCcn.ha-2D" on chromosome 2D that explains 23% of the variance. This QTL colocalized with a previously identified Cre3 locus. Novel QTL, QCcn.ha-2A, detected in the present study could be the possible unreported homeoloci to QCcn.ha-2D, QCcn.ha-2B.1 and QCcn.ha-2B.2. Six significant digenic epistatic interactions were also observed. In addition, 26 candidate genes were also identified including genes known for their involvement in PPNs (plant parasitic nematodes) resistance in different plant species. In-silico expression of putative candidate genes showed differential expression in roots during specific developmental stages. Results obtained in the present study are useful for wheat breeding to generate resistant genetic resources against H. avenae.


Phenotypic evaluation of doubled haploid lines (DHs) for H. avenae infection. In total, 114 DHs
were screened for H. avenae infection for two years. The frequency distribution for the cysts extracted is presented in Fig. 1. The frequency distribution of cyst count displayed skewed distribution towards resistance. Parental genotypes, in terms of the nematode infection in both the experiments displayed expected resistance (W-7984 or M6) and susceptibility (Opata M85) indicating that the screening was quite robust across the experiments (Supplementary Table S1). A wide range of phenotypic variability (mean ± standard deviation value of 10.56 ± 8.09) was observed among the doubled haploid (DH) lines for cysts count that ranged from 2.52 to 40.54 (CV = 0.77) ( Fig. 2; Supplementary Tables S1, S2a). There was also evidence of transgressive segregation in the population. Three DH lines displayed CCN counts (< 3.19) that were significantly (p < 0.05) lower than the resistant parent and two lines displayed susceptibility higher than the susceptible parent (> 36.40) (p < 0.05). These observations indicate the distribution or presence of favorable alleles among both parents. The heritability was estimated to be 0.83 (Supplementary Table S2b). QTL analysis. QTL analysis was performed using a high-density linkage map containing > 12 K markers (12,093 markers). Marker distribution across the wheat genome, as well as the sub-genomes, is provided in Sup-  Table 1) using CIM. These QTLs were distributed on 18 of the 21 chromosomes (except 3D, 4D and 6D). The LOD values of detected QTLs ranged from 3.02 to 11.37, with additive effects from 0.025 to 3.925 and average phenotypic variance from 0.0007 to 22.84% (Table 1). A major QTL designated as QCcn.ha-2D was detected on chromosome 2D in both years and CY (combined data) ( Table 1; Fig. 3; Supplementary Figs. S3, S4, S5, S6). This QTL was linked to marker Excalibur_rep_ c67599_2154 with LOD values ranging between 9.23 and 11.37 in the analysis of Y1, Y2 and CY and located at the distal end of chromosome 2D (136.40 centimorgan (cM)). This QTL explained 23% of the phenotypic variance in the combined year analysis with 3.925 additive effects ( Table 1). The resistance allele for this QTL was derived from the parent M6.
Another significant QTL, QCcn.ha-2A detected on chromosome 2A with LOD values ranged in between 5.30 and 9.05 at 113.40 cM in Y1, Y2 and CY analysis. This QTL was linked to marker Kukri_rep_c68068_95. This QTL explained only 3% of the total phenotypic variance in the combined year analysis with 1.628 additive effects. On chromosome 5B, another significant QTL designated as QCcn.ha-5B linked to SNP marker IACX6116 at 134. 40  www.nature.com/scientificreports/   (Table 3). No hit was found for Kukri_c29142_1046 marker located on chromosome 4A. The maximum numbers of candidate genes (six each) were detected on chromosomes 2D and 3A linked with Excalibur_rep_c67599_2154 and BS00086051_51 markers. According to the Wheat Expression Browser database, the majority of these genes were found to be expressed in wheat tissues including roots ( Fig. 6; Supplementary Tables S5, S6). Interestingly, we identified a number of genes showing differential expression in roots. We also selected nine high confidence CGs, with the majority showing more than 3 TPM (transcripts per -million) expressions in wheat roots including genes like homeobox domain, glycosyltransferase family 92 and family 8 and galacturonosyltransferase known for their involvement in PPN resistance.

Discussion
The ultimate goal of understanding the genetics of plant disease resistance is to uncover resistance loci and subsequently use this information to develop high yielding resistant cultivars. An understanding of the genetics governing resistance to CCN in wheat and the identification of molecular markers associated with the resistance gene(s) will facilitate the efforts of introgression resistance genes/QTLs into wheat. Although only 114 DH lines were used in the present study but dense genotyping allowed identification of 19 QTLs. Six of these 19 QTLs colocalized or nearby to previously reported Cre genes, QTLs or MTAs for H. avenae or H. filipjevi, while the remaining 13 QTLs were novel. Nine QTLs were detected in all three data sets. This also suggests the robustness of phenotyping and QTL analysis can be affected by several factors that include type and size of mapping, number of markers, genetic background and mapping methods 38 . As suggested in previous studies also, inconsistent QTL detection across three data sets, in the present study, might happen due to no expression or weak expression of the QTLs 39-41 .
In the present study, a highly significant QTL (QCcn.ha-2D), linked with SNP marker Excalibur_rep_ c67599_2154 was mapped on chromosome 2D. This QTL lies at 136.0-136.40 cM interval in the distal end of long arm of chromosome 2D. The resistance in D genome against H. avenae was first reported in Ae. tauschii and synthesized allo-hexaploid wheat by Eastwood et al. 15 . Several other studies reported the presence of a resistance gene against H. avenae on chromosome 2D 16,[42][43][44][45] . Interestingly, QCcn.ha-2D (detected in the present study) is co-localized with the Cre3 resistance gene (Fig. 4), which was originally identified in Ae. tauschii and introduced into cultivated wheat via a synthetic hexaploid 15,44 . By constructing two genetic linkage maps of wheat chromosome 2DL, the authors were able to develop a diagnostic microsatellite marker Xgwm301 that mapped very close to the Cre3 gene. Al-Doss et al. 46 also confirmed the presence of Cre3 gene on 2DL using Xgwm301 www.nature.com/scientificreports/ marker in another population. Eastwood et al. 16 and Ogbonnaya et al. 19 also developed an RFLP marker csE20-2 and a dominant PCR marker Cre3spf/r with complete linkage to Cre3 gene. These markers were used to identify CCN resistance in wheat cultivars 47,48 . Two QTLs identified on the distal region of long arm of 2B chromosome in the present study, namely QCcn. ha-2B.1 and QCcn.ha-2B.2 are linked to markers wsnp_RFL_Contig3917_4326857, 109.7-111.3 cM interval and Tdurum_contig12159_468, 87.50-130.0 cM interval, respectively. Both these QTLs are co-localized with the previously reported Cre1 gene (Fig. 4) and impart resistance to H. avenae in wheat 24,49,50 . This was the first gene reported for H. avenae resistance in wheat. Given that both these QTLs were discovered on chromosome 2B, they are almost certainly homeoloci of QCcn.ha-2D as also supported by a previous study 43 Fig. 4). This hypothesis is further supported by the detection of a significant MTA in the same region of the 2A chromosome (unpublished GWAS results). Further, epistatic interaction (additive × additive effect) was also observed (discussed later) between QCcn.ha-2D region and other loci on 2A (119.8 cM). However, a more detailed analysis, with additional plant material and marker enrichment for this particular chromosomal region is required to confirm this assumption. Based on the present study, it www.nature.com/scientificreports/ is evident that when the Indian pathotype was used for inoculation, detection of this novel QTL was possible, suggesting pathotype or race specific resistance expression. On chromosome 1B another minor QTL linked to Ku_c1932_1583 was identified; the physical location of this QTL is near to the QTL, QCre.srd-1B as also reported by William et al. 51 and Jayatilake et al. 52 . On chromosome 5B, significant QTL QCcn.ha-5B with LOD score ranging between 4.39 and 5.86 was identified. These QTLs co-localized with the previously identified MTA reported by Dababat et al. 53 . On chromosome 5D another minor QTL linked to wms0174 was also identified. These QTLs co-localized with the previously identified MTA reported by Mulki et al. 54 (Table 1; Fig. 4).
In the present study, the following five novel QTLs were identified: QCcn.ha-1D, QCcn.ha-3B, QCcn.ha-6B, QCcn.ha-7B and QCcn.ha-7D; these were located on chromosomes 1D, 3B, 6B, 7B and 7D. Through association mapping studies, Dababat et al. 53 and Mulki et al. 54 reported the presence of MTAs for H. avenae resistance, with different genetic and physical positions, on chromosomes 1D, 6B and 7D in synthetic hexaploid wheat and CIMMYT advanced spring wheat lines (Table 1; Fig. 4). William et al. 51,55 reported the resistance gene, Cre8 and a major QTL on the long arm of chromosome 6B against H. avenae in the DH population of Trident/Molineux. William et al. 55 also validated the linkage of H. avenae resistance QTL/gene with the RFLP marker Xcdo347-6B in another DH population of Barunga/Suneca. Jayatilke et al. 52 improved Trident/Molineux linkage map by adding 600 markers and confirmed the Cre8 locus near the distal region of chromosome 6BL as a major QTL.
Two minor QTLs, namely QCcn.ha-1A and QCcn.ha-1B, were also detected in the present study on chromosome 1A and 1B, respectively. Singh et al. 56 and William et al. 51 reported QTLs, Qcre.pau-1A and QCre.srd-1B on chromosomes 1A and 1B, respectively. In a recent study, Al-Ateeq et al. 57 a major novel resistance locus was mapped on chromosome 1BS.
Based on the published literature, 25 putative candidate genes were identified related to biotic stress or disease resistance including genes involved in PPN resistance (for details see Supplementary Table S4). For QCcn. ha-2A, two candidate genes encoding Homeobox domain, Peptidase T2 and Asparaginase2 were identified. Homeobox transcription factor is involved in maintaining the feeding sites of nematode Meloidogyne javanica infection in tomato 58 and also reported to be involved in rust resistance in wheat 59 . Two genes encoding Glycosyltransferase family 92 and family 8 were linked to QTLs QCcn.ha-7D and QCcn.ha-1A, respectively. This www.nature.com/scientificreports/ gene was reported to be involved in resistance for M. incognita (root-knot nematode) in tomato 60 , H. avenae in barley 61 and Fusarium head blight (FHB) in wheat 62 . For QCcn.ha-1A, Galacturonosyltransferase was also found. Galacturonosyltransferase is involved in susceptibility in Gossypium hirsutum at the later stage of M. incognita infection 63 . In-silico expression of identified candidate genes including homeobox domain, glycosyltransferase family 92, 8 and galacturonosyltransferase showed expression in roots at different developmental stages further suggesting a possible potential role of these candidate genes in nematode resistance. However, further delimitation of the candidate genic region is required in order to identify real candidate genes underlying the main effect QTLs in the above identified genes 64 .
Inter-locus interactions i.e., gene-by-gene interactions or epistasis and gene-by-environment interactions in QTL mapping improve the efficiency of QTL discovery. However, as also suggested in the case of soybean cyst nematode (Heterodera glycines) 65 , phenotyping for H. avenae in the present study was performed in a controlled environment and hence gene-by-environment interactions were not considered in QTL analysis. Hence, in the present study, epistatic interactions along with QTLs with additive effects were addressed for QTL mapping of H. avenae resistance. Although few studies have focused on the effect of epistasis of QTLs for soybean cyst nematode (SCN), H. glycines resistance 65-68 , root-knot nematode (RKN), M. incognita resistance 69,70 , and no epistatic interaction study has been conducted for H. avenae resistance yet. In the present study, we also conducted epistatic interaction analysis to detect the interaction of main effect QTLs with other loci (Table 2;  Supplementary Table S3). Epistatic interaction involving QCcn.ha-2D showed reduced phenotypic effect suggesting a possible application of QCcn.ha-2D alone in pre-breeding programme as a major gene which will likely show stable performance across different environments. Other main effect QTLs with smaller phenotypic effects showed improved phenotypic effect when they interact with other loci suggesting their possible utility, in combination, for MAS.
Allelic contribution may be quite effective in improving the resistance in wheat by marker-assisted breeding. Highly favorable alleles detected in the present study can be used for developing nematode resistant wheat genotypes. Alleles leading to the decrease in cyst count (Fig. 5) were considered as favorable alleles. While developing the ITMI mapping population, W7984 (M6), was used as a female parent and Opata M85 was used as a male parent 71 . W7984 (M6) and Opata M85 are originally resistant and susceptible to H. avenae pathotype infection, respectively. Interestingly, results of the present study showed that both parents contributed favourable alleles (Table 1; Fig. 6). For two stable QTLs (QCcn.ha-2D and QCcn.ha-3B), allele from the resistant parent led to the maximum decrease in cyst count (Fig. 5). Transfer of favourable alleles from the susceptible parent (Opata M85) suggests that expression of resistant favorable alleles remain masked in the susceptible parent and uncovered after recombination with the resistant parent. www.nature.com/scientificreports/

Conclusion
Identification and selection of resistant genotypes (DH lines here) with resistance loci and tolerance genetic loci may be useful for marker assisted backcrossing (MAB) to impart resistance against H. avenae. Detection of markers linked to resistance and identification of genotypes with favorable alleles for the mapped QTLs and the identified QTLs can speed up the MAB programme. The following conclusions can be drawn from the present study: except, QCcn.ha-2D, most of the significant main effect QTLs explained a small fraction of the observed phenotypic variation. Therefore, pyramiding these QTLs while taking into account positive/negative interaction may be the most effective strategy for increasing resistance to H. avenae. Eight minor QTLs (detected in all data sets) explained a combined 14% phenotypic variation. Thus, while planning a pre-breeding strategy of introducing these QTLs in pre-breeding lines, then perhaps one needs to be careful because these smaller effect QTLs may or may not be transferred. Despite having a small phenotypic effect, combining and pyramiding of these potential loci, including loci showing epistatic effects with increased phenotypic effect, for resistance into wheat breeding lines containing other sources of resistance to CCN would potentially enhance resistance level in the field 54 . Safari et al. 27 reported that Cre3 locus has the greatest impact on reducing the number of female cysts, followed by Cre1 and Cre8. Markers associated with QTLs detected in the present study, and colocalized with Cre1 and Cre3, can be used for marker assisted breeding. Development of high-throughput marker assay like KASP (competitive allele-specific PCR) followed by transfer of these QTLs to high yielding wheat lines followed by multi-location evaluation of advanced generation lines to check their agronomic performance and adaptability to diverse environments shall be the next priority. Deciphering biological function of potential candidate genes, especially those already known for their involvement in PPNs resistance/defense, is also important.

Methods
Plant material. In the present study, 114 lines, out of a total of 215 F 1 -derived doubled haploid lines (designated as SynOpDH) of the novel ITMI (International Triticeae Mapping Initiative) mapping population, were used 71 . This population was obtained from a cross between CIMMYT synthetic hexaploid wheat (SHW) line W-7984 (also known as M6) and the hard red spring wheat (HRSW) Opata M85 71 . M6 is a synthetic line www.nature.com/scientificreports/ obtained from a cross between the durum (T. turgidum L.) genotype Altar 84 and Ae. tauschii Coss., the progenitor of the bread wheat D genome 71 . Seed material of the population was originally obtained from IPK (Institute of Plant Genetics and Crop Plant Research), Gatersleben, Germany.
In the present study, the parents of the mapping population were first screened to ascertain whether this population is segregating for the resistance to H. avenae. For this purpose, we screened 10 plants of each parent in batches with the prevalent local H. avenae population from North India. Parent M6 was found to be resistant and Opata M85 was found to be susceptible. Indian wheat variety RAJ MR1 was used as a resistant check 72 and WH147, PBW343 were used as susceptible checks 73,74 in all the screening experiments. The DH lines were previously genotyped with the wheat 90 K chip and the map was prepared 71 . In total, 12,093 markers (11,678 SNPs and 415 SSRs markers) were used, after filtering the genotypic data, for the QTL mapping.

Nematode inoculums preparation and plant infection test.
Cysts were extracted from the nematode-infested soil by decanting and sieving method 75 . Cysts were surface sterilized with 0.5% sodium hypochlorite (NaOCl or bleach) for 10 min and rinsed six times with distilled water. These cysts were stored at 4 °C for two months. To enhance the hatching, J2 cysts were transferred to room temperature (~ 25 °C) for 24 h before plant inoculation.
The DH lines were screened for resistance against H. avenae under controlled environment conditions (22 °C ± 2, 16 h light, and 8 h darkness and ~ 65% relative humidity) in the growth room. For this purpose, DH lines were grown in a completely randomized design (CRD) with a minimum of five replicates in batches (over two years). The plants were irrigated with Hoagland media at a fixed interval to maintain their growth and health. The soil used was sieved and double sterilized to ensure that no other uncontrolled infection affects the plants. Seeds were pre-germinated on a wet filter paper and then transferred into 150 cm 3 PVC (polyvinyl chloride) pipes filled with steam-sterilized soil. Three holes were made in the soil around the stem base and juveniles were inoculated near the roots using a micropipette and roots were again covered with soil. Plants were inoculated with 1000 nematodes (J2) 10 days after germination and irrigated regularly.

Cyst extraction and categorization of infected plants.
White and brown cysts were extracted from both soil and roots after 75 days of infection by decanting and sieving method. Roots containing soil samples were collected in a beaker and then filled with water. The suspension was stirred well and left for about 30 s to allow the heavy sand and soil debris to settle down. Eventually, the content was poured through 850, 250 and 150 µm sieves. This process was repeated three times to ensure all cysts were collected. Also, roots were washed gently on the upper sieve to remove cysts attached to the roots. Cysts from both roots and soil were retained on a 250 µm sieve. Roots were also examined under the microscope to confirm the removal of cysts. Cyst counting was performed for each plant under a stereomicroscope (Nikon SMZ645).

Statistical analysis and QTL mapping.
In total, 2736 data points (variety, batches or environments and replication combinations) were obtained on the screened 114 lines. All data points were manually checked for any discrepancy before analysis. Later outlier Z-test was subsequently employed. After removing the outliers from the data mixed model analysis REML (restricted maximum likelihood) implemented in ASREML-R library in R was used to obtain BLUP values of the screened DH-lines. All terms are fitted as random and the best model was selected based on the change in log-likelihood differences. Finally, the best model was fitted to compute BLUPs of the DH lines that were subsequently used for QTL mapping.
QTL analysis was conducted on three data sets, namely Year 1 (Y1), Year 2 (Y2) and combined data (CY) for the two years using interval mapping. R/QTL 76 was used for interval mapping in two steps. Initially, scan one function and the three methods ("em", "hk" and "imp") implemented in R/QTL were used. All three methods detected identical QTL peaks over the chromosomes, therefore composite interval mapping (CIM) was performed using Haley-Knott regression (hk) method in R/QTL. For CIM, the number of markers as a covariate was set at three, with a window size of 5 cM. The LOD threshold value was empirically estimated with a significance level of p = 0.05. Further, we focused the study on the significant QTLs that surpassed the LOD value > 3.0. Besides, ICIM-EPI functions of the inclusive composite interval mapping (ICIM) program of QTL IciMapping v4.2 77 were employed to detect possible digenic epistatic interactions between marker loci.
Candidate gene identification and in-silico gene expression. For the identification of putative candidate genes and their gene ontology, the sequence of flanking or linked markers of the identified QTLs were blasted against the T. aestivum genome sequence information hosted at Ensembl plant release version 52 database (http:// plants. ensem bl.org/Triticum_aestivum). Important genes within the interval with known functions for disease resistance were shortlisted and considered as potential candidate genes for CCN resistance. The protein products of genes present in the flanking sequence available for the marker with maximum bases (50,000 bases before and after the marker) were reported as putative genes.
The identified putative candidate genes (CGs) were subjected to in-silico expression analysis considering different tissue and developmental stages of wheat via the 'Wheat Expression Browser-expVIP' (expression Visualization and Integration Platform; http:// www. wheat-expre ssion. com) 78,79 . Further, among all tissues, roots were selected to check the expression of candidate genes at the specific developmental stage.